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Abstract 

In this paper, a fractional generalization of the wave equation that 
describes propagation of damped waves is considered. In contrast to 
the fractional diffusion-wave equation, the fractional wave equation 
contains fractional derivatives of the same order a, 1 < a < 2 both 
in space and in time. We show that this feature is a decisive factor 
for inheriting some crucial characteristics of the wave equation like a 
constant propagation velocity of both the maximum of its fundamental 
solution and its gravity and mass centers. Moreover, the first, the 
second, and the Smith centrovelocities of the damped waves described 
by the fractional wave equation are constant and depend just on the 
equation order a. The fundamental solution of the fractional wave 
equation is determined and shown to be a spatial probability density 
function evolving in time that possesses finite moments up to the order 
a. To illustrate analytical findings, results of numerical calculations 
and numerous plots are presented. 
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1 Introduction 

During the last few decades, fractional order differential equations have been 
successfully employed for modeling of many different processes and systems, 
see e.g. |27j for different applications of derivatives and integrals of fractional 
order in physics, chemistry, engineering, astrophysics, etc. and [1] for appli- 
cations of fractional differential equations in classical mechanics, quantum 
mechanics, nuclear physics, hadron spectroscopy, and quantum field theory. 
For other interesting models in form of fractional differential equations we 
refer the reader to [3], [5], [9]- [12], [E]-[T7], [20] to mention only few of many 
recent publications. 



Among many other applications, models for anomalous transport pro- 
cesses in form of time- and/or space- fractional advection-diffusion-wave 
equations enjoyed a particular attention and have been considered by a 
number of researches since 1980's. In particular, this kind of phenomena is 
known to occur in viscoelastic media that combine characteristics of solid- 
like materials that exhibit waves propagation and fluid-like materials that 
support diffusion processes (see e.g. the recent book [32]). It is impor- 
tant to note that anomalous transport models are usually first formulated 
in stochastic form in terms of the so called continuous time random walk 
processes. The time- and/or space- fractional differential equations are then 
derived from the stochastic models for a special choice of the probability 
density functions with the infinite first or/and second moments (see e.g. 



It is well known that whereas diffusion equation describes a process, 
where a disturbance of the initial conditions spreads infinitely fast, the 
propagation velocity of the disturbance is constant for the wave equa- 
tion. In a certain sense, the time-fractional diffusion-wave equation of order 
a, 1 < a < 2 interpolates between these two different behaviors: its re- 
sponse to a localized disturbance spreads infinitely fast, but the maximum 
of its fundamental solution disperses with a finite velocity v(t,a) that is 
determined by the formula (see e.g. [6] or fH] ) 



For a = 1 (diffusion), the propagation velocity is equal to zero because of 
C\ = 0, for a = 2 (wave propagation) it remains constant and is equal to 
Ci = 1, whereas for all intermediate values of a the propagation velocity 
of the maximum point depends on time t and is a decreasing function that 
varies from +oo at time t = 0+ to zero as t — > +oo. This fact makes it 
difficult to interpret solutions to the fractional diffusion-wave equation as 
waves in case 1 < a < 2. 

In this paper, a fractional wave equation that contains fractional deriva- 
tives of the same order a, 1 < a < 2 both in space and in time is considered. 
The fractional derivative in time is interpreted in the Caputo sense whereas 
the space-fractional derivative is taken in form of the inverse operator to 
the Riesz potential (Riesz fractional derivative). It turns out that this fea- 
ture of the fractional wave equation is a decisive factor for inheriting some 
crucial characteristics of the wave equation. In particular, we show that six 
different velocities of the damped waves that are described by the funda- 
mental solution of the fractional wave equation (propagation velocity of its 



ma, [m,or [2oj). 



v(t, a) = C a t 2 
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maximum, its gravity and mass centers, the first, the second, and the Smith 
centrovelocities) are constant and depend just on the equation order a. 

From the mathematical viewpoint the fractional wave equation we deal 
with in this paper has been considered in all probability for the first time in 
[7], where an explicit formula for the fundamental solution of this equation 
was derived. In [18J, a space-time fractional diffusion- wave equation with 
the Riesz-Feller derivative of order a £ (0, 2] and skewness 9 has been inves- 
tigated in detail. A particular case of this equation called neutral-fractional 
diffusion equation that for 9 = corresponds to our fractional wave equa- 
tion has been shortly mentioned in [18]. Still, to the best knowledge of the 
author, both in-depth mathematical treatment and physical interpretation 
of the fractional wave equation seem to be not yet given in the literature. 

The rest of the paper is organized as follows. In the 2nd section, basic 
definitions, problem formulation, and some analytical results for an initial- 
value problem for a model one-dimensional fractional wave equation are 
presented. The fundamental solution G a for this problem is derived in terms 
of elementary functions for all values of a, 1 < a < 2. Moreover, G a is 
interpreted as a spatial probability density function evolving in time whose 
moments up to order a are finite. For the fundamental solution G a both its 
maximum location and its maximum value are determined in closed form. 
Remarkably, the product of the maximum location and the maximum value 
of G a is time-independent and just a function of a. Finally we show that 
both the maximum and the gravity and mass centers of the fundamental 
solution G a propagate with constant velocities like in the case of the wave 
equation but in contrast to the wave equation (a = 2) these velocities are 
different each to other for a fixed value of a, 1 < a < 2. Moreover, the first, 
the second, and the Smith centrovelocities of the damped waves described by 
the fractional wave equation are shown to be constants that depend just on 
the equation order a. To illustrate analytical findings, results of numerical 
calculations, numerous plots, their physical interpretation and discussion 
are presented in 3rd section. The last section contains some conclusions and 
open problems for further research. 
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2 Analysis of the fractional wave equation 



2.1 Problem formulation 

In this paper, we deal with the model one-dimensional fractional wave equa- 
tion 

Dfu = qT^, xeM, teM + ,l<a<2. (2) 

In (2), u = u(x,t) is a real field variable, gj^p is the Riesz space-fractional 
derivative of order a that is defined below, and Df is the Caputo time- 
fractional derivative of order a: 

(D a f)(t) := (I n - a f^)(t), n - 1< a < n, n G N (3) 

I a , a > being the Riemann-Liouville fractional integral 

{I a fm := f rR Io(t - rr-'fir) dr, a > 0, 
\f(t), a = 0, 

and r the Euler gamma function. For a = n, n £ JN, the Caputo fractional 
derivative coincides with the standard derivative of order n. 

All quantities in ^ are supposed to be dimensionless, so that the coef- 
ficient by the Riesz space-fractional derivative can be taken to be equal to 
one without loss of generality. 

For the equation ([2]), the initial- value problem 

du 

u(x,0) = <p(x) , — M) = 0, x £ M (4) 

is considered. In this paper, we are mostly interested in behavior and prop- 
erties of the fundamental solution (Green function) G a of the equation ([2]), 
i.e. its solution with the initial condition (f(x) = 5(x), 5 being the Dirac 
delta function. 

For a sufficiently well-behaved function / the Riesz space-fractional 
derivative of order a, < a < 2 is defined as a pseudo-differential operator 
with the symbol — (|23j): 

(T^f)( K ):=-\ K \ a f(K), (5) 

J- being the Fourier transform of a function /. For < a < 2, a / 1, §5§ 
can be written in the form ( |24j ) 



d\x\ a y ' 2T{-a) cos(a?r) J £ 
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For a = 1, the relation ([5]) can be interpreted in terms of the Hilbert trans- 
form 



where the integral is understood in the sense of the Cauchy principal value. 
In particular, the equation ^ with a = 1 that we call modified advection 
equation is written in the form 



that is of course different from the standard advection equation. 

For a = 2, equation ^ is reduced to the one-dimensional wave equation. 
In what follows, we focus on the case 1 < a < 2 because the case a = 2 
(wave equation) is well studied in the literature. 

2.2 Fundamental solution of the fractional wave equation 

We start our analysis by applying the Fourier transform with respect to 
the space variable x to the equation ^ with 1 < a < 2 and to the initial 
conditions Q with y?(x) = 5(x). Using definition of the Riesz fractional 
derivative, for the Fourier transform G a we get the initial-value problem 






(8) 



for the fractional differential equation 



(D a G a ){t) + \K\ a G a (K,t) = 0. 



(9) 



The unique solution of Q, Q is given by the expression (see e.g. [T5] ) 

G a (K,t) = E a {-\ K \ a t a ) 



(10) 



in terms of the Mittag-Leffler function 



OG 



k 




(11) 



The well known formula (see e.g. |21j ) 




), m £ IN, x — > +oo 
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for asymptotics of the Mittag-Leffler function that is valid for < a < 2 



and the formula (10) show that G a belongs to Li(M) with respect to k for 
1 < a < 2. Therefore we can apply the inverse Fourier transform and get 
the representation 

G a (x,t) = — / e~ iKX Ea{-\K\ a t^)dK, x e M, t > (12) 

2?r J-oo 

for the Green function G a . The last formula shows that the fundamental 
solution G a is an even function in x 

G a (-x,t) = G a (x,t), x e M, t>0 

and (12) can be rewritten as the cos- Fourier transform: 



1 f'°° 

G a [x,t) = ~ I cos(Kx)Es(-K a tP)dK, x G M, t > 0. (13) 

7T Jo 

Remarkably, the fundamental solution G a can be represented in terms of 
elementary functions for every a, 1 < a < 2. Indeed, for x = the integral 



in the right-hand side of ( 13 ) is reduced to the Mellin integral transform of 
the Mittag-Leffler function at the point s = - (for definition and proper- 
ties of the Mellin integral transform see e.g. |19j). It converges under the 
conditions a > 1 and its value is given by the formula (|19j) 

1 t'°° 1 t'°° i 

- / E a {-K a t a )dK = / E a {-u)u«~ l du (14) 

7r J Q nat J 

i r(i)rfi-i) 

v a I V a. ' 



nat rfl-a 1 



0, t > 



a ' 



because the gamma function has a pole at the point z = 0: 1/T(0) = 0. 
Since G a is an even function, we consider the integral in the right-hand 



side of (13) just in the case x = \x\ > and recognize that this integral 
is the Mellin convolution of the functions = E^{—E, a t a ) and /(£) = 
^|COs(l/^) in point 1/x. Using the known Mellin integral transforms of 
the cos-function and the Mittag-Leffler function as well as some elementary 
properties of the Mellin integral transform (|19j) we get the formulas: 

r°° 1 r (-) r fi — -) 

/(.)= I = ^ 0<RW<«, 

• f- < s » = 7^ £ ftlf' 0<s<s)<L 
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These formulas together with the convolution rule and the inverse Mellin 
integral transform lead to the representation 



G a (x, t) 



i l r* 00 r (f) r & - 1) r a - 1 ) 

7 _ ioo r(i-a) 2-r(§) 



7rax 27ri 



(is, < 7 < 1 
(15) 



of the fundamental solution G a in terms of the general Fox H-function (see, 
for example 



the form 



G a (x,t) 



1 1 

ax 2-ni 



Sj). The representation (15) can be simplified to 

7+ioo pfA^pQ _ $ 



7— JOO 



r(i-§)r(f) 



and then to the form 



G a (x,t) 



1 1 

ax 27ri 



7+ioo 



7—200 



sin(7rs/2) 
sin(7rs/a) 



ds, < 7 < a (16) 



(is, < 7 < a (17) 



by using the duplication and reflection formulas for the gamma function. 
A useful representation 



1 



G a (x,t) = -L a (t/x), x = \x\ > 0, t > 



(18) 



of the fundamental solution G a in terms of an auxiliary function L a that 
depends on the quotient t/x can be obtained from (16) or (17). Moreover, 
because (17) is in form of an inverse Mellin transform we deduce from (17) 
the Mellin transform of the function L n as follows: 



L, 



[TIT 



1 sin(7rs/2) 
a sin(7rs/a) 



(19) 



It follows from (17) that (19) holds true at least for < 3?(s) < a but as we 
see later, in fact (19) is valid even for —a < 3?(s) < a. 

Now let us represent the special case (16) of the H-function in form of 
some convergent series that can be summated in explicit form in terms of 
some elementary functions. General theory of the Mellin-Barnes integrals 
presented e.g. in [19J says that the integral in (16) is convergent under the 
condition < a < 2. For < t < x, the contour of integration in the 
integral (16) can be transformed to the loop £_oo starting and ending at 
infinity and encircling all poles s^ = — ak, k = 0,1,2, ... of the function 
r(s/a). Taking into account the relation 



res.. 



.TYs) 



(-1)* 
kl 



k = 0,1,2, 
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the residue theorem provides us with the desired series representation: 

G ^- ax 2^ kl r(-ffe)r(i-f fc ) W ( 0) 



fc=0 

that can be transformed to the form 



1 00 / t a \ k 

G a (x, t) = V sin( a7 rfe/2) - — (21) 

7TX \ X a I 

1 — 1 \ / 



fc=l 

by using the reflection formula for the gamma function. 
Now we use the summation formula 



OO / OO \ * j a s 

Vr fc S in(fa) = 3 VrV k = 9 ( ) = , (22) 

^ K ' \ ^ ) \ 1 - re m J 1 - 2r cos a + r 2 v ; 

k=l \k=l / v 7 

that is valid for a S M, \r\ < 1 to summate the series in (21 ) and obtain the 
nice representation 

c (xt)= 1 ^^H 2 ) f23l 

QV ' ; vrt 2a + 2x«i Q cos(7ra/2) +x 2a V ; 

for the Green function G a that is valid for < t < x. 

In the case < x < t we can transform the contour of integration in 
(16) to the loop £+00 encircling all poles = a(l + k), k = 0,1,2,... 
of the function r (l — — ). Applying the residue theorem we arrive at the 
representation 

1 ™ a (—l\k r(l + ifc) /j\a(i+l) 

Ga{x ' t) = ^ ho k] r(f(fc + i))r(i-f(fc + i)) \i) (24) 

that can be transformed to the form 



G a (x,t) = - — £sm(a7rfc/2) (-— ) (25) 



by using the reflection formula for the gamma function. The formula (22) 
applied to the series from the right-hand side of (25) again leads to the 
representation (23), this time for < x < t. Finally, validity of the formula 
( 23 ) for < x = t follows from the principle of analytic continuation for the 
Mellin-Barnes integrals. Thus the fundamental solution G Q for the fractional 
wave equation is given by the formula 

„ , , 1 |x| Q - 1 t a sin(vra/2) _ . . 

GJx,t) = — = -L-L / 7~r , ,„ , t > 0, x £ M (26) 

y ' 7rt 2a + 2\x\ a t a cos(ira/2) + \x\ 2a y ' 

for 1 < a < 2. 



S 



2.3 Fundamental solution as a pdf 



We begin by a remark that the formula (26) is valid for a = 1 (modified 
advection equation Q), too, that can be proved by direct calculations. In 
this case we get the well known Cauchy kernel 



Gi(x,t) 



t 



1 



7T t 2 + X 2 



(27) 



that is a spatial probability density function evolving in time. 

For a = 2 (wave equation), the Green function G2 is known to be given 
by the formula 

G 2 (x,t) = l(5(x- t) + 6(x + t)). (28) 



The representation (26) thus leads to an interesting relation 



la— 1+a 



lim 



t a sin(7ra/2) 



7T 



2-0 t 2a + 2|x| a t a cos(vra/2) + \x\ 2a 2 
for the Dirac <5-function. 



(d(x-t) + 5(x + t)), t > 0, x e M 



For 1 < a < 2, the Green function (26) is a spatial probability density 
function evolving in time, too. Indeed, the function (26) is evidently non- 
negative for all t > 0. Furthermore, for all t > and 1 < a < 2 the 
integral 



G a (x, t) dx 



2 f + °° sin(7ra/2) 
7ra Jo 1 + 2-u cos(-7ra/2) + u- 



du 



(29) 



is identically equal to 1 that can be checked by direct calculations. Thus G a 
given by ( 26 ) is a spatial probability density function evolving in time that 
can be considered to be a fractional generalization of the Cauchy kernel (27) 
for the case of an arbitrary index a, 1 < a < 2. 

Now let us study the properties of the fundamental solution (26) as a 
pdf. Because G a is an even function we can restrict our attention to the 
case x > and consider the function 



G+(x,t) 



1 



x 



a—l+a 



t a sin(7ra/2) 



7r t 2a + 2x a t a cos(7ra/2) + x 



2a 



x > 0, t > 0, 1 < a < 2. 



It is easy to see that behaves like a power function in x both at x 
and at x = +00 for a fixed t > 0: 



~-a-l „ 



0. 



+00. 



(30) 
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This means that the pdf G a possesses all finite moments up to the order a. 
In particular, the mean value of G a (its first moment) exists for all a > 1 
(we note that the Cauchy kernel does not possess a mean value). Let us 
now evaluate the moments of the one-sided fractional Cauchy kernel for 



a fixed t > 0. To do this, we refer to the representation (18) of in terms 



of the auxiliary function L a that can be now represented in the form 
_ , . 1 T°sin(7ra/2) 

L a (t ) = = y ' -, r > 0, 1< a < 2. (31 

y; TTT 2a + 2r Q cos(vra/2) + 1' v ; 

Taking into account this formula, the function prL a (|T|) can be interpreted 
as a fractional Cauchy pdf of order a. 

The moment of the order /3, \f3\ < a of G^ can be represented in terms 



of the Mellin integral transform of L a that is known (see the formula (19)) 
and thus evaluated: 

rGUx,t)xPdx = te f Mrlr^dr^ ^f) . (32) 
Jo Jo a sin(7r/3/a) 

In particular, we get the formula 

1 

G+{x,t)dx = - (33) 



o 



that is in accordance with (29) because G a is an even function. 



We mention also important formula 



oc 







G±(x,t)xdx = t ——-,l<a<2 (34) 

asm(7r/a) 



for the first moment of the one-sided fractional Cauchy kernel . 

2.4 Extrema points, gravity and mass centers of G a , and 
location of its energy 

Now we derive some important analytical properties of the fractional Cauchy 



kernel (26). First we remark that G a (0,t) = and G a (x,t) > for x ^ 0, 
so that x = is a minimum point for the fundamental solution G a for any 
t > 0. Because G a is an even function we again consider its restriction to 
x > 0, i.e. the function G+. 

To determine locations of maxima of G^ for fixed values of t and a we 
solve the equation 
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that turns out to be equivalent to the quadratic equation 



(a + 1) 



with solutions given by 

x a 



+ 2cos(7ra/2) 



(a - 1) = 



cos(7ra/2) ± \J a 2 - sin 2 (7ra/2) 
a + 1 



Since we are interested in nonnegative solutions the only candidate for this 
role is the point 



x 



cos(7ra/2) + \J a 2 — sin 2 (7ra/2) 



a + 1 



(35) 



Because ^^(x,t) is positive for < c a and negative for > c a we 
conclude that the point 

Xa(t) = v p (a)t, v p (a) := (c a )s (36) 



with c a given by (35) is the only maximum point of the fractional Cauchy 
kernel G a for x > 0. Of course, the point — x*(t) < is another maximum 
point of G a because G a is an even function. 

To determine the maximum value of the function G a that coincides with 
the maximum value of and is denoted by G* (t) we substitute the point 
x = x* (t) given by ( 36 ) into the function and get 



i 



c Q sin(7ra/2) 



TTV p (a)t 1 + 2c a cos(7ra/2) + c 2 ' 



(37) 



where v p (a) and c a are defined as in the formulas (35) and (36). 

It follows from the formulas (36) and (37) that for a fixed value of a, 1 < 
a < 2 the product p a of the maximum value G* a {t) and the maximum 
location x* a (t) is time-independent: 



Po 



c a sin(7ra/2) 



7T 1 + 2c a cos(7ra/2) + c 2 



(38) 



This means that the maximum point (x* (t), G* (t)) of the fundamental so- 
lution G a with a fixed value of a, 1 < a < 2 moves in time along a hyperbola 
that is completely determined by the value of a (see the formula (38)). It 
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is interesting to note that the product p a coincides with the value of the 
fundamental solution G a at the point (l,v p (a)): 



Pa 



1 c a sin(7ra;/2) 

7r 1 + 2c a cos(7ra/2) + c a 



G a (l,v p (a)). 



This property can be probably used to give a physical interpretation of the 
formula (38). 

Now we calculate the location of the gravity center x a (t) of the funda- 
mental solution G a that is defined by the formula (we recall that G a is an 
even function) 

Jq 00 x G a (x, t) dx 



J °° G a (x,t) dx 



Using the formulas (33) and ( |34[ ) we get the following result: 

a sm{ir/a) 

The "mass" -center x™(t) of G a is determined by the formula ([8]) 

f™ xGl(x,t)dx 



(39) 



(40) 



f™G 2 a {x,t)dx ' 



(41) 



Substituting (26) into (41) and transforming the obtained integrals we get 
the representation 

x a (*) = v m (a) t, v m (a) 



POO 

Jo T 



L 2 Jr)dT 



^T- 2 Ll{ T )dr 



(42) 



where the function L a is defined by ( |31[ ). Because the Mellin transform of 
L a is known (see the formula (19)), the integrals 



POO 

/ t 13 L\(t) dr, -2a - 1< < 2a - 1 
J o 



(43) 



can be interpreted as Mellin convolutions of the functions r /3+1 L Q ,(r) and 
L a (l/r) in point x = 1 and thus expressed in terms of an H- function with 
parameters depending on a and (3 and with the argument x = 1. Because 
there are no routines for numerical evaluation of the H-function available, 
we prefer to stay by the representation of v m (a) given in (42) and not to 
transform it to a quotient of two H- functions. Remarkably, there exists an 
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explicit formula for the integrals ( 43 ) in the case a = 1 (modified advection 
equation), namely ([22]) 

^ t?L\(t) dr = 1 + /? 1 -3</3<l. (44) 

Jo 4?r cos(vr/3/2) 

The "mass" -center x™ of Gi can be then represented by the simple formula 

x?(t) = -t. (45) 

Finally we mention that "location of energy" of the damped wave that 
is represented by the fundamental solution G a is given by the formula ([2J) 

Jo Gi(x,t)dt 
and can be represented in the form 

v c (a) J TL 2 a (r)dT 

where the function L a is defined by (31 ). Because the integral J °° r L^(t) dr 
diverges for a = 1, v c {a) tends to as a -4 1. 

2.5 The velocities of the damped waves 

It is well known (see e.g. pQ, [2J, [8], or [26]) that several different definitions 
of the wave velocities and in particular of light velocity can be introduced. 
For the damped waves that are described by the fractional wave equation 
([2]) we calculate propagation velocity of the maximum of its fundamental so- 
lution G a that can be interpreted as the phase velocity, propagation velocity 
of the gravity center of G a , the velocity of its "mass" -center or the pulse 
velocity, and three different kinds of its centrovelocity. It turns out that all 
these velocities are constant in time and depend just on the order a of the 
fractional wave equation. Whereas four out of six velocities are different 
each to other, the fist centrovelocity coincides with the Smith centrovelocity 
and the the second centrovelocity is the same as the pulse velocity. 

We start with the phase velocity and determine it using the formula 



(36) that leads to the result that the maximum location of the fundamental 
solution G a propagates with a constant velocity v p (a) that is given by the 
expression 



, . dxt(t) ( - cos(7ra/2) + J a 2 - sin 2 (7ra/2) \ a 

V ^ a):= ^r = [ • (48) 
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For a = 1 (modified advection equation y}), the propagation velocity of 
the maximum of G a is equal to zero (the maximum point stays at x = 0) 
whereas for a = 2 (wave equation) the maximum point propagates with the 
constant velocity 1. 

To determine the propagation velocity v g (a) of the gravity center of G a 



we employ the formula (40) and get the following result: 



v g (a) :-- 



dXa(t) 

dt 



a sin(7r/a) 



(49) 



v g (a) is thus time-independent and determined by the order a of the frac- 
tional wave equation. Evidently, v g (2) = 1 and v g {a) — > +oo as a — > 1 + 0. 

The velocity v m {a) of the "mass" -center of G a or its pulse velocity 
is obtained from the formula (|42|) and is equal to 



[a 




(50) 



the pulse velocity is 



where the function L a is defined b y <\3 1j 
equal to - ~ 0.64 (see the formula M5I) 

Following [2] we define the second centrovelocity V2{pt) as the mean pulse 



velocity computed from to time t . It follows from ( 42 ) and ( 50 ) that for the 



damped wave that is described by the fundamental solution of the fractional 
wave equation the second centrovelocity is equal to its pulse velocity v m (a): 



v 2 {a 



l (t) 



t 



JO T L a 



(r) dr 



Jo T ~ 



2 Ll(r)dr 



(51) 



The Smith centrovelocity v c {a) ([26]) of the damped waves describes the 
motion of the first moment of their energy distribution and can be evaluated 
in explicit form using the formula (14 71): 



v c (u) 



dt* 



dx 



-i 



I[Ll(r)dr ^ 

j™TL 2 a {T)dT 



(52) 



where the function L a is defined by (31 ). Because the integral r -Z^(t) dr 
diverges for a = 1, the Smith centrovelocity tends to as a — > 1. 

Finally, we calculate the first centrovelocity v\{a) that is defined as the 



mean centrovelocity from to x ( [2] ) . It follows from ( 47 ) and ( 52 ) that for 
the damped wave G a the first centrovelocity is equal to the Smith centrov- 
elocity v c (a): 



vi(a) :-- 



t%(x) 



( v_ f~L*(T)dT 

~ C{ ) ~ J °°tLUt) dr 



(53) 
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Figure 1: Fundamental solution G a : Plots for values of a = 1.01 (1st line, 
left), 1.1 (1st line, right), 1.5 (2nd line left), and 1.9 (2nd line, right) for 
-0.5 < x < 0.5 and t = 0.1, 0.2, 0.3 



As we have seen, all velocities introduced above are constant in time 
and depend just on the order a of the fractional wave equation. The phase 
velocity, the velocity of the gravity center of G a , the pulse velocity, and the 
Smith centrovelocity are different each to other whereas the fist centroveloc- 
ity coincides with the Smith centrovelocity and the second centrovelocity is 
the same as the pulse velocity. For the physical interpretation and meaning 
of the velocities that were determined above we refer to e.g. PQ, [2], or jS]. 



3 Discussion of the obtained results and plots 

To start with, let us consider evolution of the fundamental solution G a 
in time for some characteristic values of a. In Fig. [T] plots of G a for 
a = 1.01, 1.1, 1.5, and 1.9 are presented. As we can see, in all cases maxi- 



mum location is moved linearly in time according to the formula ( 36 ) whereas 



the maximum value decreases according to the formula (37). The behavior 
of G a can be thus interpreted as propagation of damped waves whose am- 
plitude decreases with time. This phenomena can be very clearly recognized 
on the 3D plot presented in Fig. [2] Of course, because of the nonlocal char- 
acter of the fractional derivatives in the fractional wave equation solutions 
to this equation show some properties of diffusion processes, too. In par- 
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0.5 



Figure 2: Plots of G a for a = 1.5, -0.5 < x < 0.5, and < t < 0.3 (left) 
and of the maximum value m a of the function G a at the time instant t = 1 
(right) 



ticular, the fundamental solution G a is positive for all x ^ at any small 
time instance t > that means that a disturbance of the initial conditions 
spreads infinitely fast and the equation ([2]) is non relativistic like the clas- 
sical diffusion equation. But in contrast to the diffusion equation, both the 
maximum of the fundamental solution G a , its gravity and mass centers and 
location of its energy propagate with the finite constant velocities like in the 
case of the fundamental solution of the wave equation. 

Plots of the propagation velocity v p of the maximum of the fundamen- 
tal solution G a (phase velocity), the velocity v g of its gravity center, its 
pulse velocity v m and its centrovelocity v c are presented in Fig. [3| As ex- 
pected, v p = v c = and v m « 0.64 for a = 1 (modified advection equation) 
and all velocities smoothly approach the value 1 as a — > 2 (wave equa- 
tion). For 1 < a < 2, v p , v m , and v c monotonously increase whereas v g 
monotonously decreases. It is interesting to note that for all velocities the 
property dv }^ (2 — 0) = holds true, i.e. in a small neighborhood of the point 
a = 2 the velocities of G a are nearly the same as in the case of the fundamen- 
tal solution of the wave equation. The velocity v g of the gravity center of G a 
tends to +oo for a. — )■ 1 + and t > (modified advection equation) because 
the first moment of the Cauchy kernel (27) does not exist. It is interesting 
to note that for all a, 1 < a < 2 the velocities v p , v g , v m , v c are different 
each to other and fulfill the inequalities v c (a) < v p (a) < v m (a) < v g (a). 
For a = 2, all velocities are equal to 1. 

In Fig. [2] (right), a plot of the maximum value of the fundamental solu- 
tion G a that is given by the formula ( 37 ) is presented for t = 1 . Surprisingly, 
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Figure 3: Plots of the phase velocity v p (a) and the gravity center velocity 
v g (a) (left) and of the pulse velocity v m (a) and the centrovelocity v c (a) 
(right) 



the function m a = G* a (1) is not monotone. Numerical calculations show that 
it has a minimum that is located at the point a ~ 1.13. At a = 1, G a (x, 1) 
has a (local) maximum with the value - ~ 0.32. Then G a (x, 1) monotoni- 



cally decreases to m r 



0.28 at a r 



1.13 (minimum) and then starts 



to monotonically increase and tends to +oo when a — > 2 — that is in ac- 
cordance with behavior of the fundamental solution (^(S(x — t) + 5(x + t)) 
of the wave equation (a = 2). We note that the location x* a (t) of the maxi- 
mum value G* a (i) of G a is given by the formula ( 36 ) . For t = 1 , location of 
the maximum value coincides with the phase velocity v p that is presented 
in Fig. [3J 

Finally, some plots of the product p a of the maximum value G* (t) and 



the maximum location x* (t) given by the formula (38) are presented in Fig. 
[4j The function p a = p a {ot) is monotone increasing for 1 < a < 2 and varies 
from at the point a = 1 (the location of the maximum of the Cauchy kernel 
( |27| ) is always at the point x = 0) to +oo at the point a = 2 (the maximum 
value of the the Green function (^(5(x — t) + 5(x + t)) of the wave equation is 
+oo). In this sense, the product p a can be considered to be a characteristic 
property of the damped waves that are described by the fractional wave 
equation Q. As we can see on the left plot, the product p a varies between 
and 10 for 1 < a < 1.98 and changes very slowly for 1.01 < a < 1.9. For 
a — > 1 + (see the right plot of Fig. [4j and a — > 2 — 0, p a goes to and +oo, 
respectively, very fast. For a damped wave that is supposed to be described 
by the fractional wave equation ([2]) with an unknown exponent a, the value 
of the product p a can be measured and used to recover a either from the 
formula (38) or graphically from Fig. |4j 



17 




1 1.2 1.4 1.6 1.8 2 1 1.01 1.02 1.03 1.04 1.05 



Figure 4: Plots of the product p a of the maximum value G* a {t) and the 
maximum location x* (t) on the interval 1 < a < 2 (left) and on the interval 
1 < a < 1.05 (right) 



4 Conclusions and open questions 

In this paper, a model one-dimensional fractional wave equation with the 
fractional derivatives of order a, 1 < a < 2 both in space and in time is 
introduced and considered in detail. The fractional wave equation inherits 
some crucial characteristics of the wave equation like constant propagation 
velocities of the maximum of its Green function, its gravity and mass centers, 
and its energy location. Because the maximum value of the fundamental 
solution G a (wave amplitude) decreases with time whereas its location moves 
with a constant velocity, solutions to the fractional wave equation can be 
interpreted as damped waves. Moreover, G a that turns out to be expressed 
in terms of elementary functions for all values of a, 1 < a < 2 can be 
interpreted as a spatial probability density function evolving in time whose 
moments up to order a are finite. For the fundamental solution G a , both its 
maximum location and its maximum value are determined in closed form. 
Remarkably, the product of the maximum location and the maximum value 
of G a is time-independent and just a function of a. 

Among problems for further research we mention two- and three- 
dimensional fractional wave equations with different initial or / and boundary 
conditions. Of course, it would be interesting to consider fractional wave 
equations with fractional derivatives defined in different ways. We note here 
that in [18] a space-time fractional diffusion-wave equation with the Caputo 
derivative of order (5 G (0, 2] in time and the Riesz- Feller derivative of order 
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a 6 (0, 2] and skewness 9 in space has been investigated in detail. A par- 
ticular case of this equation called neutral-fractional diffusion equation that 
for 9 = corresponds to our fractional diffusion equation has been shortly 
mentioned in |18j . 

Another interesting and important problem for further research would be 
determination of other velocities like the group velocity or the ratio-of-units 
velocity (see e.g. [I] or [8]) for the damped waves described by the fractional 
wave equations and comparison them each to other at least for the linear 
equations with the constant coefficients. Finally, fractional wave equations 
with non-constant coefficients as well as qualitative behavior of solutions of 
non-linear fractional wave equations would be worth to consider. 

Acknowledgment: The author is thankful to Prof. Francesco Mainardi 
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during author's visit to the University of Bologna in December 2011. 
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